#analysis for the appendix 

#Appendix A: codebook: no code------------------


#Appendix B: Sample Characteristics---------------------------
library(dplyr)
library(psych)
library(tidyr)

#for experiment 1
selected_columns <- data_correct[, c("RACE", "GENDER", "PID", "IDEOLOGY.", "REL_IMPORTANCE", "INCOME", "AGE", "CN.Score")]

sapply(selected_columns, class)

#edit selected columns as needed as numeric
selected_columns <- selected_columns %>% 
  mutate(RACE = case_when(
    RACE == "White, non- Hispanic"  ~ 1, 
    RACE == "Hispanic " ~ 2, 
    RACE == "Asian" ~ 3, 
    RACE == "Black or African American, non- Hispanic" ~ 4, 
    RACE == "American Indian or Alaskan Native" ~ 5, 
    RACE ==  "Multiple races, non-Hispanic" ~ 6, 
    RACE == "Other:" ~ 7, 
    RACE == "Native Hawaiian or other Pacific Islander" ~ 8
  ))

selected_columns <- selected_columns %>%
  mutate(GENDER = case_when( 
    GENDER == "Man" ~ 1, 
    GENDER == "Non-binary / third gender-" ~ 2, 
    GENDER == "Prefer not to say" ~ 3, 
    GENDER == "Woman" ~ 4))


selected_columns <- selected_columns %>%
  mutate(INCOME = case_when( 
    INCOME == "Less than $10,000" ~ 1, 
    INCOME == "$10,000 - $19,999" ~ 2, 
    INCOME == "$20,000 - $29,999" ~ 3, 
    INCOME == "$30,000 - $39,999" ~ 4, 
    INCOME ==  "$40,000 - $49,999" ~ 5, 
    INCOME == "$50,000 - $59,999" ~ 6, 
    INCOME == "$60,000 - $69,999" ~ 7, 
    INCOME ==  "$70,000 - $79,999" ~ 8, 
    INCOME ==  "$80,000 - $89,999" ~ 9, 
    INCOME ==  "$90,000 - $99,999" ~ 10,
    INCOME == "$100,000 - $149,999" ~ 11, 
    INCOME == "More than $150,000" ~12))

selected_columns <- selected_columns %>%
  mutate(AGE = case_when( 
    AGE == "18 - 24" ~ 1, 
    AGE == "25 - 34" ~ 2, 
    AGE == "35 - 44" ~ 3, 
    AGE == "45 - 54" ~ 4, 
    AGE ==  "55 - 64" ~ 5, 
    AGE == "65 - 74" ~ 6, 
    AGE == "75 - 84" ~ 7))

#run CN.Score as numeric

# Descriptive statistics using dplyr
summary_stats <- selected_columns %>%
  summarise(
    mean_RACE = mean(RACE, na.rm = TRUE),
    sd_RACE = sd(RACE, na.rm = TRUE),
    mean_GENDER = mean(GENDER, na.rm = TRUE),
    sd_GENDER = sd(GENDER, na.rm = TRUE),
    mean_PID = mean(PID, na.rm = TRUE),
    sd_PID = sd(PID, na.rm = TRUE),
    mean_IDEOLOGY. = mean(IDEOLOGY., na.rm = TRUE),
    sd_IDEOLOGY. = sd(IDEOLOGY., na.rm = TRUE),
    mean_REL_IMPORTANCE = mean(REL_IMPORTANCE, na.rm = TRUE),
    sd_REL_IMPORTANCE = sd(REL_IMPORTANCE, na.rm = TRUE),
    mean_INCOME = mean(INCOME, na.rm = TRUE),
    sd_INCOME = sd(INCOME, na.rm = TRUE),
    mean_AGE = mean(AGE, na.rm = TRUE),
    sd_AGE = sd(AGE, na.rm = TRUE),
    mean_CN.Score = mean(as.numeric(CN.Score), na.rm = TRUE),
    sd_CN.Score = sd(as.numeric(CN.Score), na.rm = TRUE)
  )

summary_stats

#counts

# Convert all relevant columns to character to resolve data type mismatch
data_correct <- data_correct %>%
  mutate(across(c(RACE, RELIGION, GENDER, PID, IDEOLOGY., 
                  REL_IMPORTANCE, AGE, block_id, CN.Score),
                as.character))

# Define a lookup table for variable names and their labels
variable_labels <- tibble(
  Variable = c("RACE", "RELIGION", "GENDER", "PID", "IDEOLOGY.", 
               "REL_IMPORTANCE", "AGE", "block_id", "CN.Score"),
  Description = c("Race", "Religion", "Gender", "PID", "Ideology", 
                  "Importance of Religion", "Age", 
                  "CN Measurement Bin", "CN Score")
)

# Pivot the data into a long format, count occurrences, and join with labels
counts <- data_correct %>%
  pivot_longer(cols = c(RACE, RELIGION, GENDER, PID, IDEOLOGY., 
                        REL_IMPORTANCE, AGE, block_id, CN.Score),
               names_to = "Variable", values_to = "Value") %>%
  count(Variable, Value, name = "Count") %>%
  left_join(variable_labels, by = "Variable")

# Arrange for easier viewing
counts <- counts %>%
  arrange(Variable, desc(Count))


# View the result
print(counts)

#experiment 1 durability check 
selected_columns_dur <- dur_data_correct[, c("RACE", "GENDER", "PID", "IDEOLOGY.", "REL_IMPORTANCE", "INCOME", "AGE", "CN.Score")]

sapply(selected_columns_dur, class)

#edit selected columns as needed as numeric
selected_columns_dur <- selected_columns_dur %>% 
  mutate(RACE = case_when(
    RACE == "White, non- Hispanic"  ~ 1, 
    RACE == "Hispanic " ~ 2, 
    RACE == "Asian" ~ 3, 
    RACE == "Black or African American, non- Hispanic" ~ 4, 
    RACE == "American Indian or Alaskan Native" ~ 5, 
    RACE ==  "Multiple races, non-Hispanic" ~ 6, 
    RACE == "Other:" ~ 7, 
    RACE == "Native Hawaiian or other Pacific Islander" ~ 8
  ))

selected_columns_dur <- selected_columns_dur %>%
  mutate(GENDER = case_when( 
    GENDER == "Man" ~ 1, 
    GENDER == "Non-binary / third gender-" ~ 2, 
    GENDER == "Prefer not to say" ~ 3, 
    GENDER == "Woman" ~ 4))


selected_columns_dur <- selected_columns_dur %>%
  mutate(INCOME = case_when( 
    INCOME == "Less than $10,000" ~ 1, 
    INCOME == "$10,000 - $19,999" ~ 2, 
    INCOME == "$20,000 - $29,999" ~ 3, 
    INCOME == "$30,000 - $39,999" ~ 4, 
    INCOME ==  "$40,000 - $49,999" ~ 5, 
    INCOME == "$50,000 - $59,999" ~ 6, 
    INCOME == "$60,000 - $69,999" ~ 7, 
    INCOME ==  "$70,000 - $79,999" ~ 8, 
    INCOME ==  "$80,000 - $89,999" ~ 9, 
    INCOME ==  "$90,000 - $99,999" ~ 10,
    INCOME == "$100,000 - $149,999" ~ 11, 
    INCOME == "More than $150,000" ~12))

selected_columns_dur <- selected_columns_dur %>%
  mutate(AGE = case_when( 
    AGE == "18 - 24" ~ 1, 
    AGE == "25 - 34" ~ 2, 
    AGE == "35 - 44" ~ 3, 
    AGE == "45 - 54" ~ 4, 
    AGE ==  "55 - 64" ~ 5, 
    AGE == "65 - 74" ~ 6, 
    AGE == "75 - 84" ~ 7))

#run CN.Score as numeric

# Descriptive statistics 
summary_stats_dur <- selected_columns_dur %>%
  summarise(
    mean_RACE = mean(RACE, na.rm = TRUE),
    sd_RACE = sd(RACE, na.rm = TRUE),
    mean_GENDER = mean(GENDER, na.rm = TRUE),
    sd_GENDER = sd(GENDER, na.rm = TRUE),
    mean_PID = mean(PID, na.rm = TRUE),
    sd_PID = sd(PID, na.rm = TRUE),
    mean_IDEOLOGY. = mean(IDEOLOGY., na.rm = TRUE),
    sd_IDEOLOGY. = sd(IDEOLOGY., na.rm = TRUE),
    mean_REL_IMPORTANCE = mean(REL_IMPORTANCE, na.rm = TRUE),
    sd_REL_IMPORTANCE = sd(REL_IMPORTANCE, na.rm = TRUE),
    mean_INCOME = mean(INCOME, na.rm = TRUE),
    sd_INCOME = sd(INCOME, na.rm = TRUE),
    mean_AGE = mean(AGE, na.rm = TRUE),
    sd_AGE = sd(AGE, na.rm = TRUE),
    mean_CN.Score = mean(as.numeric(CN.Score), na.rm = TRUE),
    sd_CN.Score = sd(as.numeric(CN.Score), na.rm = TRUE)
  )

summary_stats_dur

# Convert all relevant columns to character to resolve data type mismatch
dur_data_correct <- dur_data_correct %>%
  mutate(across(c(RACE, RELIGION, GENDER, PID, IDEOLOGY., 
                  REL_IMPORTANCE, AGE, block_id, CN.Score),
                as.character))

# Define a lookup table for variable names and their labels
variable_labels_dur <- tibble(
  Variable = c("RACE", "RELIGION", "GENDER", "PID", "IDEOLOGY.", 
               "REL_IMPORTANCE", "AGE", "block_id", "CN.Score"),
  Description = c("Race", "Religion", "Gender", "PID", "Ideology", 
                  "Importance of Religion", "Age", 
                  "CN Measurement Bin", "CN Score")
)

# Pivot the data into a long format, count occurrences, and join with labels
counts_dur <- dur_data_correct %>%
  pivot_longer(cols = c(RACE, RELIGION, GENDER, PID, IDEOLOGY., 
                        REL_IMPORTANCE, AGE, block_id, CN.Score),
               names_to = "Variable", values_to = "Value") %>%
  count(Variable, Value, name = "Count") %>%
  left_join(variable_labels, by = "Variable")

# Arrange for easier viewing
counts_dur <- counts_dur %>%
  arrange(Variable, desc(Count))

# View the result
print(counts_dur)


#experiment 2
selected_columns_2 <- data2[, c("race", "gender", "pid", "ideology", 
                                "rel_importance", "income", "age", "CN.Score2")]

sapply(selected_columns_2, class)

#edit selected columns as needed as numeric
selected_columns_2 <- selected_columns_2 %>% 
  mutate(race = case_when(
    race == "White, non- Hispanic"  ~ 1, 
    race == "Hispanic " ~ 2, 
    race == "Asian" ~ 3, 
    race == "Black or African American, non- Hispanic" ~ 4, 
    race == "American Indian or Alaskan Native" ~ 5, 
    race ==  "Multiple races, non-Hispanic" ~ 6, 
    race == "Other:" ~ 7, 
    race == "Native Hawaiian or other Pacific Islander" ~ 8
  ))

selected_columns_2 <- selected_columns_2 %>%
  mutate(gender = case_when( 
    gender == "Man" ~ 1, 
    gender == "Non-binary / third gender-" ~ 2, 
    gender == "Prefer not to say" ~ 3, 
    gender == "Woman" ~ 4))


selected_columns_2 <- selected_columns_2 %>%
  mutate(income = case_when( 
    income == "Less than $10,000" ~ 1, 
    income == "$10,000 - $19,999" ~ 2, 
    income == "$20,000 - $29,999" ~ 3, 
    income == "$30,000 - $39,999" ~ 4, 
    income ==  "$40,000 - $49,999" ~ 5, 
    income == "$50,000 - $59,999" ~ 6, 
    income == "$60,000 - $69,999" ~ 7, 
    income ==  "$70,000 - $79,999" ~ 8, 
    income ==  "$80,000 - $89,999" ~ 9, 
    income ==  "$90,000 - $99,999" ~ 10,
    income == "$100,000 - $149,999" ~ 11, 
    income == "More than $150,000" ~12))

selected_columns_2 <- selected_columns_2 %>%
  mutate(age = case_when( 
    age == "18 - 24" ~ 1, 
    age == "25 - 34" ~ 2, 
    age == "35 - 44" ~ 3, 
    age == "45 - 54" ~ 4, 
    age ==  "55 - 64" ~ 5, 
    age == "65 - 74" ~ 6, 
    age == "75 - 84" ~ 7))

#run CN.Score as numeric

# Descriptive statistics
summary_stats_2 <- selected_columns_2 %>%
  summarise(
    mean_race = mean(race, na.rm = TRUE),
    sd_race = sd(race, na.rm = TRUE),
    mean_gender = mean(gender, na.rm = TRUE),
    sd_gender = sd(gender, na.rm = TRUE),
    mean_pid = mean(pid, na.rm = TRUE),
    sd_pid = sd(pid, na.rm = TRUE),
    mean_ideology = mean(ideology, na.rm = TRUE),
    sd_ideology = sd(ideology, na.rm = TRUE),
    mean_rel_importance = mean(rel_importance, na.rm = TRUE),
    sd_rel_importance = sd(rel_importance, na.rm = TRUE),
    mean_income = mean(income, na.rm = TRUE),
    sd_income = sd(income, na.rm = TRUE),
    mean_age = mean(age, na.rm = TRUE),
    sd_age = sd(age, na.rm = TRUE),
  )

summary_stats_2

# Convert all relevant columns to character to resolve data type mismatch
data2 <- data2 %>%
  mutate(across(c(race, religion, gender, pid, ideology, 
                  rel_importance, age, CN.Score2),
                as.character))

# Define a lookup table for variable names and their labels
variable_labels_2 <- tibble(
  Variable = c("race", "religion", "gender", "pid", "ideology", 
               "rel_importance", "age", "CN.Score2"),
  Description = c("Race", "Religion", "Gender", "PID", "Ideology", 
                  "Importance of Religion", "Age", "block ID")
)

# Pivot the data into a long format, count occurrences, and join with labels
counts_2 <- data2 %>%
  pivot_longer(cols = c(race, religion, gender, pid, ideology, 
                        rel_importance, age, CN.Score2),
               names_to = "Variable", values_to = "Value") %>%
  count(Variable, Value, name = "Count") %>%
  left_join(variable_labels, by = "Variable")

# Arrange for easier viewing
counts_2 <- counts_2 %>%
  arrange(Variable, desc(Count))

# View the result
print(counts_2)


#experiment 2 durability check 

selected_columns_2_dur <- data2_dur[, c("race", "gender", "pid", "ideology", 
                                        "rel_importance", "income", "age", "CN.Score2")]

sapply(selected_columns_2_dur, class)

#edit selected columns as needed as numeric
selected_columns_2_dur <- selected_columns_2_dur %>% 
  mutate(race = case_when(
    race == "White, non- Hispanic"  ~ 1, 
    race == "Hispanic " ~ 2, 
    race == "Asian" ~ 3, 
    race == "Black or African American, non- Hispanic" ~ 4, 
    race == "American Indian or Alaskan Native" ~ 5, 
    race ==  "Multiple races, non-Hispanic" ~ 6, 
    race == "Other:" ~ 7, 
    race == "Native Hawaiian or other Pacific Islander" ~ 8
  ))

selected_columns_2_dur <- selected_columns_2_dur %>%
  mutate(gender = case_when( 
    gender == "Man" ~ 1, 
    gender == "Non-binary / third gender-" ~ 2, 
    gender == "Prefer not to say" ~ 3, 
    gender == "Woman" ~ 4))


selected_columns_2_dur <- selected_columns_2_dur %>%
  mutate(income = case_when( 
    income == "Less than $10,000" ~ 1, 
    income == "$10,000 - $19,999" ~ 2, 
    income == "$20,000 - $29,999" ~ 3, 
    income == "$30,000 - $39,999" ~ 4, 
    income ==  "$40,000 - $49,999" ~ 5, 
    income == "$50,000 - $59,999" ~ 6, 
    income == "$60,000 - $69,999" ~ 7, 
    income ==  "$70,000 - $79,999" ~ 8, 
    income ==  "$80,000 - $89,999" ~ 9, 
    income ==  "$90,000 - $99,999" ~ 10,
    income == "$100,000 - $149,999" ~ 11, 
    income == "More than $150,000" ~12))

selected_columns_2_dur <- selected_columns_2_dur %>%
  mutate(age = case_when( 
    age == "18 - 24" ~ 1, 
    age == "25 - 34" ~ 2, 
    age == "35 - 44" ~ 3, 
    age == "45 - 54" ~ 4, 
    age ==  "55 - 64" ~ 5, 
    age == "65 - 74" ~ 6, 
    age == "75 - 84" ~ 7))

#run CN.Score as numeric

# Descriptive statistics using dplyr
summary_stats_2_dur <- selected_columns_2_dur %>%
  summarise(
    mean_race = mean(race, na.rm = TRUE),
    sd_race = sd(race, na.rm = TRUE),
    mean_gender = mean(gender, na.rm = TRUE),
    sd_gender = sd(gender, na.rm = TRUE),
    mean_pid = mean(pid, na.rm = TRUE),
    sd_pid = sd(pid, na.rm = TRUE),
    mean_ideology = mean(ideology, na.rm = TRUE),
    sd_ideology = sd(ideology, na.rm = TRUE),
    mean_rel_importance = mean(rel_importance, na.rm = TRUE),
    sd_rel_importance = sd(rel_importance, na.rm = TRUE),
    mean_income = mean(income, na.rm = TRUE),
    sd_income = sd(income, na.rm = TRUE),
    mean_age = mean(age, na.rm = TRUE),
    sd_age = sd(age, na.rm = TRUE),
  )

summary_stats_2_dur

# Convert all relevant columns to character to resolve data type mismatch
data2_dur <- data2_dur %>%
  mutate(across(c(race, religion, gender, pid, ideology, 
                  rel_importance, age, CN.Score2),
                as.character))

# Define a lookup table for variable names and their labels
variable_labels_2_dur <- tibble(
  Variable = c("race", "religion", "gender", "pid", "ideology", 
               "rel_importance", "age", "CN.Score2"),
  Description = c("Race", "Religion", "Gender", "PID", "Ideology", 
                  "Importance of Religion", "Age", 
                   "block ID")
)

# Pivot the data into a long format, count occurrences, and join with labels
counts_2_dur <- data2_dur %>%
  pivot_longer(cols = c(race, religion, gender, pid, ideology, 
                        rel_importance, age, CN.Score2),
               names_to = "Variable", values_to = "Value") %>%
  count(Variable, Value, name = "Count") %>%
  left_join(variable_labels, by = "Variable")

# Arrange for easier viewing
counts_2_dur <- counts_2_dur %>%
  arrange(Variable, desc(Count))

# View the result
print(counts_2_dur)

#Appendix C: Full Hypotheses: no code ------------------------
#Appendix D: Results from Full Sample------------------------

#match what's in the paper: no controls, interaction with continuous variable 

full_model <- estimatr::lm_robust(CN_IDENTIFICATION ~ Treatment + as.numeric(CN.Score), 
                                  data = data)
summary(full_model)

#Public
#Creating scale for public identification: using questions FRIEND_EVENT_CN, CN_WE, FURTHER_CONVERSATION
data$public_id_scale <-  rowMeans(data[, c("FRIEND_EVENT_CN", "CN_WE", "FUTHER_CONVERSATION")], 
                                  na.rm = TRUE)


full_model_public <- estimatr::lm_robust(public_id_scale ~ Treatment + as.numeric(CN.Score), 
                                         data = data)
summary(full_model_public)



#output tables
full_list <- list(full_model)
texreg(full_list, include.ci = FALSE)


full_list_public <- list(full_model_public)
texreg(full_list_public, include.ci = FALSE)



#Appendix E: Removal of Control for White Nationalism: no code -----------------------

#Appendix F: Original Results from Durability Check--------------------
#packages: 
library(dplyr)
library(ggplot2)



#Subsetting into treatment blocks


#create data frames for CN blocks 
dur_data_rejecters_correct <- dur_data_correct[dur_data_correct$block_id == 1, ]
dur_data_skeptics_correct <- dur_data_correct[dur_data_correct$block_id == 2, ]
dur_data_sympathizers_correct <- dur_data_correct[dur_data_correct$block_id == 3, ]
dur_data_adherents_correct <- dur_data_correct[dur_data_correct$block_id== 4, ]



#Regressions with controls for durability check--------------------

# Run the linear regression model
dur_model <- estimatr::lm_robust(DURABILITY_3 ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                 data = dur_data_correct, fixed_effects = ~block_id)
summary(dur_model)
linhyp_dur <- car::linearHypothesis(dur_model, "TreatmentP = TreatmentN")
linhyp_dur


#rejecters
dur_model_rejecters <- estimatr::lm_robust(DURABILITY_3 ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                           data = dur_data_rejecters_correct, fixed_effects = ~block_id)
summary(dur_model_rejecters)
#comparison of P and N
linhyp_dur_rejecters <- car::linearHypothesis(dur_model_rejecters, "TreatmentP = TreatmentN")
linhyp_dur_rejecters

#skeptics
dur_model_skeptics <- estimatr::lm_robust(DURABILITY_3 ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                          data = dur_data_skeptics_correct, fixed_effects = ~block_id)
summary(dur_model_skeptics)

#comparison of P and N
linhyp_dur_skeptics <- car::linearHypothesis(dur_model_skeptics, "TreatmentP = TreatmentN")
linhyp_dur_skeptics


#sympathizers
dur_model_sympathizers <- estimatr::lm_robust(DURABILITY_3 ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                              data = dur_data_sympathizers_correct, fixed_effects = ~block_id)
summary(dur_model_sympathizers)

#comparison of P and N
linhyp_dur_sympathizers <- car::linearHypothesis(dur_model_sympathizers, "TreatmentP = TreatmentN")
linhyp_dur_sympathizers

#adherents
dur_model_adherents <- estimatr::lm_robust(DURABILITY_3 ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                           data = dur_data_adherents_correct, fixed_effects = ~block_id)
summary(dur_model_adherents)

#comparison of P and N
linhyp_dur_adherents <- car::linearHypothesis(dur_model_adherents, "TreatmentP = TreatmentN")
linhyp_dur_adherents





#Public results with controls for durability check-------------------
#Creating scale for public identification: using questions FRIEND_EVENT_CN, CN_WE, FURTHER_CONVERSATION
dur_data_correct$public_id_scale <-  rowMeans(dur_data_correct[, c("DURABILITY_1.", "DURABILITY_2", "DURABILITY_4")], na.rm = TRUE)
dur_data_rejecters_correct$public_id_scale<-  rowMeans(dur_data_rejecters_correct[, c("DURABILITY_1.", "DURABILITY_2", "DURABILITY_4")], na.rm = TRUE)
dur_data_skeptics_correct$public_id_scale<-  rowMeans(dur_data_skeptics_correct[, c("DURABILITY_1.", "DURABILITY_2", "DURABILITY_4")], na.rm = TRUE)
dur_data_sympathizers_correct$public_id_scale<-  rowMeans(dur_data_sympathizers_correct[, c("DURABILITY_1.", "DURABILITY_2", "DURABILITY_4")], na.rm = TRUE)
dur_data_adherents_correct$public_id_scale<-  rowMeans(dur_data_adherents_correct[, c("DURABILITY_1.", "DURABILITY_2", "DURABILITY_4")], na.rm = TRUE)




dur_model_public <- estimatr::lm_robust(public_id_scale ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                        data = dur_data_correct, fixed_effects = ~block_id)
summary(dur_model_public)

#comparison of P and N
linhyp_dur_public <- car::linearHypothesis(dur_model_public, "TreatmentP = TreatmentN")
linhyp_dur_public

#rejecters
dur_model_rejecters_public <- estimatr::lm_robust(public_id_scale ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                                  data = dur_data_rejecters_correct, fixed_effects = ~block_id)
summary(dur_model_rejecters_public)
#comparison of P and N
linhyp_dur_rejecters_public <- car::linearHypothesis(dur_model_rejecters_public, "TreatmentP = TreatmentN")
linhyp_dur_rejecters_public

#skeptics
dur_model_skeptics_public <- estimatr::lm_robust(public_id_scale ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                                 data = dur_data_skeptics_correct, fixed_effects = ~block_id)
summary(dur_model_skeptics_public)

#comparison of P and N
linhyp_dur_skeptics_public <- car::linearHypothesis(dur_model_skeptics_public, "TreatmentP = TreatmentN")
linhyp_dur_skeptics_public

#sympathizers
dur_model_sympathizers_public <- estimatr::lm_robust(public_id_scale ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                                     data = dur_data_sympathizers_correct, fixed_effects = ~block_id)
summary(dur_model_sympathizers_public)

#comparison of P and N
linhyp_dur_sympathizers_public <- car::linearHypothesis(dur_model_sympathizers_public, "TreatmentP = TreatmentN")
linhyp_dur_sympathizers_public


#adherents

dur_model_adherents_public <- estimatr::lm_robust(public_id_scale ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                                  data = dur_data_adherents_correct, fixed_effects = ~block_id)
summary(dur_model_adherents_public)

#comparison of P and N
linhyp_dur_adherents_public <- car::linearHypothesis(dur_model_adherents_public, "TreatmentP = TreatmentN")
linhyp_dur_adherents_public


#output tables
dur_list <- list(dur_model, dur_model_adherents, dur_model_sympathizers, dur_model_skeptics, dur_model_rejecters)
texreg(dur_list, include.ci = FALSE)

dur_list_public <- list(dur_model_public, dur_model_adherents_public, dur_model_sympathizers_public, dur_model_skeptics_public, dur_model_rejecters_public)
texreg(dur_list_public, include.ci = FALSE)

#Appendix G: ENE: no code --------------


#Appendix H: Full treatment language: no code ----------


#Appendix I: Original Results as Preregistered --------------
#results as preregistered,

model <- estimatr::lm_robust(CN_IDENTIFICATION ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3,
                             fixed_effects = ~block_id, data = data_correct)
summary(model)

linhyp_model <- car::linearHypothesis(model, "TreatmentP = TreatmentN")
linhyp_model


#models of each CN measurement bin for comparison
#create data frames for CN blocks 
data_rejecters_correct <- data_correct[data_correct$block_id == 1, ]
data_skeptics_correct <- data_correct[data_correct$block_id == 2, ]
data_sympathizers_correct <- data_correct[data_correct$block_id == 3, ]
data_adherents_correct <- data_correct[data_correct$block_id == 4, ]



#Rejecters- Block ID/Bin 1

model_rejecters <- estimatr::lm_robust(CN_IDENTIFICATION ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                       data = data_rejecters_correct, fixed_effects = ~block_id)
summary(model_rejecters)

linhyp_rejecters <- car::linearHypothesis(model_rejecters, "TreatmentP = TreatmentN")
linhyp_rejecters
#Skeptics- Block ID/Bin 2
model_skeptics <- estimatr::lm_robust(CN_IDENTIFICATION ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                      data = data_skeptics_correct, fixed_effects = ~block_id)
summary(model_skeptics)


linhyp_skeptics <- car::linearHypothesis(model_skeptics, "TreatmentP = TreatmentN")
linhyp_skeptics
#Sympathizers- Block ID/Bin 3
model_sympathizers <- estimatr::lm_robust(CN_IDENTIFICATION ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                          data = data_sympathizers_correct, fixed_effects = ~block_id)
summary(model_sympathizers)

linhyp_sympathizers <- car::linearHypothesis(model_sympathizers, "TreatmentP = TreatmentN")
linhyp_sympathizers

#Adherents- Block ID/Bin 4
model_adherents <- estimatr::lm_robust(CN_IDENTIFICATION ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                       data = data_adherents_correct, fixed_effects = ~block_id)
summary(model_adherents)

linhyp_adherents <- car::linearHypothesis(model_adherents, "TreatmentP = TreatmentN")
linhyp_adherents

#output tables
# Positive vs negative - Ha1 and Ha2
main_list <- list(model, model_adherents, model_sympathizers, model_skeptics, model_rejecters)
texreg(main_list, include.ci = FALSE)

#for public identification-----

data_rejecters_correct$public_id_scale<-  rowMeans(data_rejecters_correct[, c("FRIEND_EVENT_CN", "CN_WE", "FUTHER_CONVERSATION")], na.rm = TRUE)
data_skeptics_correct$public_id_scale<-  rowMeans(data_skeptics_correct[, c("FRIEND_EVENT_CN", "CN_WE", "FUTHER_CONVERSATION")], na.rm = TRUE)
data_sympathizers_correct$public_id_scale<-  rowMeans(data_sympathizers_correct[, c("FRIEND_EVENT_CN", "CN_WE", "FUTHER_CONVERSATION")], na.rm = TRUE)
data_adherents_correct$public_id_scale<-  rowMeans(data_adherents_correct[, c("FRIEND_EVENT_CN", "CN_WE", "FUTHER_CONVERSATION")], na.rm = TRUE)


#Regressions for public measurement 

#Creating scale for public identification: using questions FRIEND_EVENT_CN, CN_WE, FURTHER_CONVERSATION
data_correct$public_id_scale <-  rowMeans(data_correct[, c("FRIEND_EVENT_CN", "CN_WE", "FUTHER_CONVERSATION")], na.rm = TRUE)
data_rejecters_correct$public_id_scale<-  rowMeans(data_rejecters_correct[, c("FRIEND_EVENT_CN", "CN_WE", "FUTHER_CONVERSATION")], na.rm = TRUE)
data_skeptics_correct$public_id_scale<-  rowMeans(data_skeptics_correct[, c("FRIEND_EVENT_CN", "CN_WE", "FUTHER_CONVERSATION")], na.rm = TRUE)
data_sympathizers_correct$public_id_scale<-  rowMeans(data_sympathizers_correct[, c("FRIEND_EVENT_CN", "CN_WE", "FUTHER_CONVERSATION")], na.rm = TRUE)
data_adherents_correct$public_id_scale<-  rowMeans(data_adherents_correct[, c("FRIEND_EVENT_CN", "CN_WE", "FUTHER_CONVERSATION")], na.rm = TRUE)

#Regressions

# Run the linear regression model
model_public <- estimatr::lm_robust(public_id_scale ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                    data = data_correct, fixed_effects = ~block_id)
summary(model_public)

linhyp_public <- car::linearHypothesis(model_public, "TreatmentP = TreatmentN")
linhyp_public

#Looking at comparison between positive and negative

#models of each CN measurement bin for comparison


#Rejecters- Block ID/Bin 1

model_rejecters_public <- estimatr::lm_robust(public_id_scale ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                              data = data_rejecters_correct, fixed_effects = ~block_id)
summary(model_rejecters_public)

linhyp_rejecters_public <- car::linearHypothesis(model_rejecters_public, "TreatmentP = TreatmentN")
linhyp_rejecters_public

#Skeptics- Block ID/Bin 2
model_skeptics_public <- estimatr::lm_robust(public_id_scale ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                             data = data_skeptics_correct, fixed_effects = ~block_id)
summary(model_skeptics_public)

linhyp_skeptics_public <- car::linearHypothesis(model_skeptics_public, "TreatmentP = TreatmentN")
linhyp_skeptics_public

#Sympathizers- Block ID/Bin 3
model_sympathizers_public <- estimatr::lm_robust(public_id_scale ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                                 data = data_sympathizers_correct, fixed_effects = ~block_id)
summary(model_sympathizers_public)

linhyp_sympathizers_public <- car::linearHypothesis(model_sympathizers_public, "TreatmentP = TreatmentN")
linhyp_sympathizers_public

#Adherents- Block ID/Bin 4
model_adherents_public <- estimatr::lm_robust(public_id_scale ~ Treatment + race_binary + PID + REL_IMPORTANCE + CN_ISSUE_3, 
                                              data = data_adherents_correct, fixed_effects = ~block_id)
summary(model_adherents_public)


linhyp_adherents_public <- car::linearHypothesis(model_adherents_public, "TreatmentP = TreatmentN")
linhyp_adherents_public

#output tables
# Positive vs negative - Ha1 and Ha2
public_list <- list(model_public, model_adherents_public, model_sympathizers_public, model_skeptics_public, model_rejecters_public)
texreg(public_list, include.ci = FALSE)

#packages: 
library(dplyr)
library(ggplot2)
library(DeclareDesign)
library(texreg)
library(cobalt)
library(WeightIt)
#data visualizations

# Function to calculate mean and confidence intervals
calc_mean_ci <- function(data_correct, conf_level = 0.95) {
  n <- length(data_correct)
  mean_val <- mean(data_correct)
  stderr <- sd(data) / sqrt(n)
  error_margin <- qt(conf_level / 2 + 0.5, n - 1) * stderr
  lower_bound <- mean_val - error_margin
  upper_bound <- mean_val + error_margin
  return(c(mean = mean_val, lower = lower_bound, upper = upper_bound))
}

# Summarize the data to calculate mean, confidence intervals, and counts
data_summary <- data_correct %>%
  filter(!is.na(Treatment)) %>%  # Remove rows with NA in Treatment
  group_by(block_id, Treatment) %>%
  summarise(mean = mean(CN_IDENTIFICATION),
            n = n(), # Count the number of observations
            ci_lower = ifelse(n > 1, mean(CN_IDENTIFICATION) - qt(0.975, df = n - 1) * sd(CN_IDENTIFICATION) / sqrt(n), NA),
            ci_upper = ifelse(n > 1, mean(CN_IDENTIFICATION) + qt(0.975, df = n - 1) * sd(CN_IDENTIFICATION) / sqrt(n), NA),
            .groups = 'drop')  # Prevents creation of a grouped data frame

# Plot with labels
ggplot(data_summary, aes(x = block_id, y = mean, fill = Treatment)) +
  geom_bar(position = position_dodge(width = 0.9), stat = "identity", alpha = 0.7) +
  geom_text(aes(label = n), 
            position = position_dodge(width = 0.9), 
            vjust = 9.5, size = 3, hjust = 0.5) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), 
                position = position_dodge(width = 0.9), 
                width = 0.25, colour = "red", alpha = 0.9) + 
  labs(x = "Christian Nationalism Scale Group", y = "Mean Identification with Christian Nationalism") +
  scale_x_discrete(breaks = 1:4, labels = c("Rejecters", "Skeptics", "Sympathizers", "Adherents")) +
  scale_fill_manual(name = "Treatment", 
                    values = c("C" = "gray20", "P" = "gray50", "N" = "gray80"),
                    labels = c("Control" = "C", "Positive" = "P", "Negative" = "N")) +
  guides(fill = guide_legend(override.aes = list(pattern = c("none", "stripe", "checkerboard")))) + # Specify patterns
  theme_minimal()



#public identification 
# Function to calculate mean and confidence intervals
calc_mean_ci_public <- function(data_correct, conf_level = 0.95) {
  n <- length(data_correct)
  mean_val <- mean(data_correct)
  stderr <- sd(data_correct) / sqrt(n)
  error_margin <- qt(conf_level / 2 + 0.5, n - 1) * stderr
  lower_bound <- mean_val - error_margin
  upper_bound <- mean_val + error_margin
  return(c(mean = mean_val, lower = lower_bound, upper = upper_bound))
}

# Summarize the data to calculate mean, confidence intervals, and counts
data_summary_public <- data_correct %>%
  filter(!is.na(Treatment)) %>%  # Remove rows with NA in Treatment
  group_by(block_id, Treatment) %>%
  summarise(mean = mean(public_id_scale),
            n = n(), # Count the number of observations
            ci_lower = ifelse(n > 1, mean(public_id_scale) - qt(0.975, df = n - 1) * sd(public_id_scale) / sqrt(n), NA),
            ci_upper = ifelse(n > 1, mean(public_id_scale) + qt(0.975, df = n - 1) * sd(public_id_scale) / sqrt(n), NA),
            .groups = 'drop')  # Prevents creation of a grouped data frame

# Plot with labels
ggplot(data_summary_public, aes(x = block_id, y = mean, fill = Treatment)) +
  geom_bar(position = position_dodge(width = 0.9), stat = "identity", alpha = 0.7) +
  geom_text(aes(label = n), 
            position = position_dodge(width = 0.9), 
            vjust = 9.5, size = 3, hjust = 0.5) +
  geom_errorbar(aes(ymin = ci_lower, ymax = ci_upper), 
                position = position_dodge(width = 0.9), 
                width = 0.25, colour = "red", alpha = 0.9) + 
  labs(x = "Christian Nationalism Scale Group", y = "Mean Public Identification with Christian Nationalism") +
  scale_x_discrete(breaks = 1:4, labels = c("Rejecters", "Skeptics", "Sympathizers", "Adherents")) +
  scale_fill_manual(name = "Treatment", 
                    values = c("C" = "gray20", "P" = "gray50", "N" = "gray80"),
                    labels = c("Control" = "C", "Positive" = "P", "Negative" = "N")) +
  guides(fill = guide_legend(override.aes = list(pattern = c("none", "stripe", "checkerboard")))) + # Specify patterns
  theme_minimal()

#Appendix J: factor analysis --------------
# Load necessary libraries
library(psych)
library(dplyr)

#factor analysis: make sure the variables all map onto each other well--------
# Load necessary libraries
library(dplyr)
library(psych)

# Calculate the public_id_scale using rowMeans
data2 <- data2 %>%
  mutate(public_id_scale_2 = rowMeans(across(c(friend_event_cn., cn_we, further_conversation), as.numeric), na.rm = TRUE))

# Perform Factor Analysis on the selected variables
fa_result_2 <- data2 %>%
  select(friend_event_cn., cn_we, further_conversation) %>%
  na.omit() %>%
  fa(nfactors = 1, rotate = "none", fm = "ml")

# View results
print(fa_result_2)

# Optionally visualize the factor loadings
fa.diagram(fa_result_2)


#what about when you remove the friend event one? 
# Perform Factor Analysis without FRIEND_EVENT_CN
fa_result_no_friend_event_2 <- data2 %>%
  select(cn_we, further_conversation) %>%
  na.omit() %>%
  fa(nfactors = 1, rotate = "none", fm = "ml")

# View the results
print(fa_result_no_friend_event_2)

# Optionally visualize the factor loadings
fa.diagram(fa_result_no_friend_event_2)

# Experiment 1: Calculate Cronbach's alpha for the selected variables
alpha_result <- data_correct %>%
  select(FRIEND_EVENT_CN, CN_WE, FUTHER_CONVERSATION) %>%
  na.omit() %>%
  psych::alpha()
print(alpha_result)


# Experiment 1: Cronbach's alpha without FRIEND_EVENT_CN
alpha_result_no_friend_event <- data_correct %>%
  select(CN_WE, FUTHER_CONVERSATION) %>%
  na.omit() %>%
  psych::alpha()
print(alpha_result_no_friend_event)

# Experiment 2: Calculate Cronbach's alpha for the selected variables
alpha_result2 <- data2 %>%
  select(friend_event_cn., cn_we, further_conversation) %>%
  na.omit() %>%
  psych::alpha()
print(alpha_result2)

# Experiment 2: Cronbach's alpha without FRIEND_EVENT_CN
alpha_result_no_friend_event2 <- data2 %>%
  select(cn_we, further_conversation) %>%
  na.omit() %>%
 psych::alpha()
print(alpha_result_no_friend_event2)

#results with other scale 
#scale creation
#without friend_event_cn
data_correct$public_id_scale_2 <-  rowMeans(data_correct[, c("CN_WE", "FUTHER_CONVERSATION")], na.rm = TRUE)
data2$public_id_scale_2 <-  rowMeans(data2[, c("cn_we", "further_conversation")], na.rm = TRUE)

#experiment 1
model_public_nc_without <- estimatr::lm_robust(public_id_scale_2 ~ Treatment+ as.numeric(CN.Score), 
                                       data = data_correct)
summary(model_public_nc_without)

#experiment 2
model_public_nc_2_without <- estimatr::lm_robust(public_id_scale_2 ~ Treatment+ as.numeric(CN.Score), 
                                         data = data2)
summary(model_public_nc_2_without)

#outputs 
public_list_nc_without <- list(model_public_nc_without, model_public_nc_2_without)
texreg(public_list_nc_without, include.ci = FALSE)



# Appendix K: distribution of CN score by treatment groups ----------------------
data2 <- data2 %>%
  filter(!is.na(Treatment) & Treatment != "")

ggplot(data2, aes(x = CN.Score, fill = Treatment)) +
  geom_histogram(alpha = 0.6, position = "dodge", bins = 30) +
  facet_wrap(~ Treatment, scales = "free_y") +
  theme_minimal() +
  labs(title = "Distribution of CN.Score by Treatment", x = "CN.Score", y = "Frequency")



#QQ plots
ggplot(data2, aes(sample = CN.Score)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~ Treatment) +
  theme_minimal() +
  labs(title = "Q-Q Plots for CN.Score by Treatment")



#Appendix M: interaction models 
model_nc_int <- estimatr::lm_robust(CN_IDENTIFICATION ~ Treatment * as.numeric(CN.Score),
                                    + data = data_correct)
model_public_nc_int <- estimatr::lm_robust(public_id_scale ~ Treatment * as.numeric(CN.Score), 
                                       data = data_correct)
model_nc_2_int <- estimatr::lm_robust(as.numeric(cn_identification) ~ Treatment * as.numeric(CN.Score),
                                  data = data2)
model_public_nc_2_int <- estimatr::lm_robust(public_id_scale2 ~ Treatment * as.numeric(CN.Score), data = data2)
                                         
int_list <- list(model_nc_int, model_public_nc_int, model_nc_2_int, model_public_nc_2_int)

texreg(int_list, include.ci=FALSE)



#Appendix N: tests of active vs public id ------------
#for experiment 1
data_correct <- data_correct %>%
  mutate(
    active_id = as.numeric(CN_IDENTIFICATION),  # active identity
    public_id = public_id_scale                 # public identity scale
  )

data_correct_long <- data_correct %>%
  select(participant_id = ID.., Treatment, CN.Score, active_id, public_id) %>%
  pivot_longer(
    cols = c(active_id, public_id),
    names_to = "ID_type",
    values_to = "ID_outcome"
  )

#run the model
model_combined1 <- lm_robust(ID_outcome ~ Treatment * ID_type + as.numeric(CN.Score),
                                 data = data_correct_long)
summary(model_combined1)


#for experiment 2
data2_long <- data2 %>%
  mutate(
    active_id = as.numeric(cn_identification),
    public_id = public_id_scale2
  ) %>%
  select(participant_id = id, Treatment, CN.Score, active_id, public_id) %>%
  pivot_longer(
    cols = c(active_id, public_id),
    names_to = "ID_type",
    values_to = "ID_outcome"
  )

#run the model
model_combined2 <- lm_robust(ID_outcome ~ Treatment * ID_type + as.numeric(CN.Score),
                            data = data2_long)
summary(model_combined2)


#output 
appendixn_list <- list(model_combined1, model_combined2)
texreg(appendixn_list, include.ci = FALSE)

